Finite difference lattice Boltzmann model with flux limiters for liquid-vapor systems 
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In this paper we apply a finite difference lattice Boltzmann model to study the phase separation 
in a two-dimensional liquid-vapor system. Spurious numerical effects in macroscopic equations 
are discussed and an appropriate numerical scheme involving flux limiter techniques is proposed 
to minimize them and guarantee a better numerical stability at very low viscosity. The phase 
separation kinetics is investigated and we find evidence of two different growth regimes depending 
on the value of the fluid viscosity as well as on the liquid-vapor ratio. 
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I. INTRODUCTION 

Lattice Boltzmann (LB) models approach physical 
phenomena in fluid systems using a phase-space dis- 
cretizcd form of the Boltzmann equation P, U |3j S 13 ■ 
Conservation equations are derived by calculating mo- 
ments of various order of this equation H, 0, 0, 0. 
EH H^ . After the publication of the first LB model 
which exhibits phase separation LB models were 

widely used to investigate the complex behavior of single- 
or multi-component /phase fluid systems |3, |5|| and refer 



mainly to isothermal systems ^ ES 113 Ealll HI El 
122 . 123L I24 . |25| . This limitation comes from the constant 
value of the lattice speed q which in LB models is related 
to the temperature T, the lattice spacing Ss and the time 
step dt through two separate relations 
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where Cs = ^JksT /m is the isothermal speed of sound 
for an ideal fluid, m is the mass of fluid particles, % is a 
constant depending on the geometry of the lattice, and 
kB is Boltzmann's constant pi l2fil| . 

According to the "collide and stream" philosophy of 
LB models, fluid collides in the lattice nodes and there- 
after moves along the lattice links in a lapse 5t towards 
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neighboring nodes with the speed ci given by Eq. (jJl 
0, ,3, 4, 5.]. Such a relationship is no longer consid- 
ered in finite difference lattice Boltzmann (FDLB) mod- 
els mUlllllMIIJl which start directly from the Boltz- 
mann equation and have a better numerical stability. In 
such models there is more freedom to choose the discrete 
velocity set, as done recently in the thermal FDLB model 
of Watari and Tsutahara [321 where the possibility of hav- 
ing different sets of velocities allows to release the con- 
straint of constant temperature. Also, the use of FDLB 
models is promising, e.g., when considering LB models 
for multicomponent fluid systems, where the masses of 
the components are not identical and Eq. would lead 
to different lattice speeds. In this context, FDLB models 
may be viewed as a convenient alternative to interpola- 
tion supplemented LB models ^3, 34, 35]. 

FDLB models, as well as LB models, are known to in- 
troduce spurious terms in the mass and momentum con- 
servation equations, which are dependent on the lattice 
spacing 5s and the time step 5t [31|. The behavior of an 
isothermal fluid system subjected to FDLB simulation is 
governed by the apparent values of the viscosity and/or 
diffusivity. The expression of these quantities with re- 
spect to 5s and 5t depends on the finite difference scheme 
used in the FDLB model. Consequently, the choice of the 
numerical scheme may alter significantly the macroscopic 
behavior of the fluid system observed during simulations 
as well as the numerical stability. This problem still lacks 
necessary clariflcation and should be always considered 
in order to recover the correct physical interpretation of 
simulation results. 

The purpose of this paper is to investigate these numer- 
ical aspects by using a FDLB model addressing the phase 
separation kinetics in a van der Waals fluid. Phase sepa- 
ration in liquid-vapor systems has not received as much 
attention as in binary fluids ,;3S] . Under the hypothesis 
of dynamical scaling the late time kinetics can be charac- 
terized in terms of a single length scale R{t) which grows 
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according to the pow er law R{t) ^ t", where a is the 
growth exponent p37| . The late time growth, when hy- 
drodynamics is neglected, is expected to be described by 
the AUcn-Cahn theory which gives a growth exponent 
a = 1/2 js^. When hydrodynamics comes into play, 
the liquid-vapor system behaves similarly to binary flu- 
ids so that a growth exponent a = 2/3 is expected [36|. 
Previous numerical studies used molecular dynamics sim- 
ulations Isaliol an d a LB model based on a free energy 
functional |4ll 143. li^ . In molecular dynamics simula- 
tions it was found evidence for the growth exponent 1/2 
[3^ |40|. Bv using the free-energy LB model, Osborn 
et al. |4j] found the growth exponents 2/3 and 1/2 at 
low and high viscosity, respectively, independently on the 
system composition. Mecke and Sofonea 42, %3|, using 
the same algorithm for an off-symmetric system, found 
a crossover from 2/3 to 1/2 at low viscosity, and 1/3 at 
high viscosity. We will compare results of our model with 
the aforementioned ones. 

To model the liquid-vapor system, a standard force 
term 0, O |25| is added to the discretized Boltzmann 
equations. The resulting FDLB model is described in 
Section ^ In Section IIIII we introduce two numeri- 
cal schemes, namely, the first order upwind finite dif- 
ference scheme and a higher order one which uses flux 
limitcrs 44, There we show the difference between 
FDLB, "collide and stream" LB and volumetric LB mod- 
els 0, 01 . We thereafter discuss the spurious numer- 
ical effects these schemes introduce in the fluid equa- 
tions. Section IIVI reports the simulation results, where 
special attention was given to the effects of the numer- 
ical schemes on estimation of the growth exponent. In 
order to clarify the phenomenology and estimate accu- 
rately the exponent a, we monitored the size of domains 
R{t) by using three independent measures. A discussion 
about the method and results ends the paper. 



II. THE MODEL 

The 2D FDLB model follows the LB model for non- 
ideal fluids m El m S m. The starting point is 
provided by the set of M partial derivatives equations re- 
sulting from the discretization of the Boltzmann equation 
on a square lattice C when the collision term is linearized 
using the BGK approximation 50]. In non-dimensional 
form, this set reads 



dtfi + eipdpfi 
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i = 0,1,...7V 



(3) 



Since we will deal with a van der Waals fluid, we used the 
following reference quantities for particle number den- 
sity, temperature and speed to get the non-dimensional 
form ^ of the discretized Boltzmann equations: = 
Na/V^c: Tr = T^, cr = ^kBTjm. Here Na is Avo- 
gadro's number, Vmc is the molar volume at the critical 



point and Tc is the critical temperature. With this choice 
of reference quantities, the dimensionless speed is [sij 

c = ci/cR = ^/ehi (4) 

where 9 = T/Tr is the dimensionless temperature and 
the constant x equals 1/3 for the square lattice we use 
(see later for details on the lattice) [31|. If we take the 
system size as the reference length Ir, we get the ref- 
erence time Ir from the condition IrCr = Ir. The 
non-dimensionalized lattice spacing is deflned by the 
number of lattice nodes N we choose along the non- 
dimensionalized system length L 
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The particle distribution functions fi = fi{x, t) are de- 
fined in the nodes x of the square lattice C. In the D2Q9 
model we use in this paper, J\f — 8 and the velocities 
areiiii 
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The equilibrium distribution functions f^"^ = f^'^{x,t) 
are expressed as series expansions of the Maxwellian dis- 
tribution function, up to second order with respect to the 
local velocity u = u{x,t), whose Cartesian components 
are Ufs (521] : 
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The weight coefficients are: 
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The local density n = n(a;, f), as well as the components 
of the local velocity u which enter Eq. , are calculated 
from the distribution functions as follows: 
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The force term in Eqs. (PJ is given by [i^, |4i 



xecm 



Fp = -d^ip'^pn + lidpiV'n) (11) 



where 



and 
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are the non-dimensionahzed pressures of the ideal and the 
van der Waals fluid, respectively [sj. With the equation 
of state in the form (|13|) . the critical point is located at 
6 = 1 and n = 1. The parameter k controls the surface 
tension psf. The mass and momentum equations are 
recovered from Eqs. ||2J) after using the standard Chap- 
man - Enskog procedure up to second order with respect 
to Knudsen number Kn = ct/L. These equations read 

— EF 



(14) 
-dc^p"^ + KndaC^^n) (15) 
vdp [n {daUp + dpUa)] 



dtn + dp{nui3) 
dt{nua) + dpinuaup) 



where 



XC T 



(16) 



is the physical value of the kinematic viscosity j3ll |. The 
particular numerical scheme used to solve Eqs. Q may 
introduce a spurious viscosity term that adds to the phys- 
ical value, as seen in the next section. Finally, we note 
that the force term Ijlll) allows to recover the Navier - 
Stokes equation H15|l where the pressure appearing on 
the r.h.s. is subjected to the van der Waals equation of 
state (^ni- 



f,{x,t + St) = Mx,t)~ (17) 
c5t 

— [fiix,t) - fi{x- Ssejc,t)] + StQi{x,t) 

OS 

(18) 



where 
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n{x, t) 



dp [fix,t)-p^{x,t)]+Kdp [V^nix,t)] 



x/r'(a=> t) hf3 - up{x, t)] - - [fiix, t) - f'r'^{x, t)] , 

T 

i = 0, 1,...7V 

As discussed in Ref. finite difference schemes in- 
troduce spurious numerical terms in the conservation 
equations. This happens because the real evolution equa- 
tions recovered (up to second order in space and time) 
from the updating rules l(T7|) are 



where 
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We get the following form of the conservation equations 
up to second order in the Knudsen number: 

dtn + dp{nup) = {t/j - (f>)dadp [xc^nSap + nUaUp] (23) 



III. FINITE DIFFERENCE SCHEMES 

A. First-order upwind scheme 

The set of phase space discretized equations may 
be solved numerically by using an appropriate finite dif- 
ference scheme defined on the lattice C Simple second- 
order schemes like the centered one or the Lax - Wendroff 
scheme 0, [s^ are unstable because of large values 
of the density gradient which may occur in the interface 
regions separating the liquid and vapor phases of the van 
der Waals fluid. The first-order upwind scheme, which is 
also used in LB models S S j is a good candidate 
because of its stability. When associated to the forward 
time stepping rule, this scheme gives the following up- 
dating rule for the distribution functions defined in node 



dt{nua) + dp{nUaUp) = —daP^ + Knda{y'^n) (24) 
+ i^apdp [n {daUp + dpUa)] 
+ xc^('0 - (l))dp [da{nup) + updau + Uadpn] 

Thus, the finite difference scheme introduces spurious 
terms, depending on the quantity [ip—cj)), in both the con- 
servation equations (compare with Eqs. H14 |l - H15|l 'l. while 
the physical value Hl()|) of the kinematic viscosity is re- 
placed by the apparent value 31] 

vap^xc'ir + i^) (25) 

One could use 5s, St, and c such that ^ = cj) 5s = cSt 
and remove spurious terms in the Eqs. (|23|) - (|24|l . In this 
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FIG. 1: Characteristics lines on the square lattice, for the 
directions ei (a) and es (b). 



case it would be Uap = xc^('''+'5t/2). In order to maintain 
the apparent value of the viscosity close to the physical 
one and allow very small values of one should require 
8t <C T. Since the condition = is equivalent to ask 
N = L/cSt, one should have iV - 10"* when r ~ lO'^, be- 
ing c L 1. This would require a huge computational 
effort when doing 2D or 3-D simulations using the first 
order upwind FDLB model. Higher order flux limiter 
schemes provide a possibility to overcome this problem 
giving a better stability. As we will see further, these 
schemes improve the accuracy of the FDLB simulations 
with respect to the upwind scheme, for the same value of 
the number N of lattice nodes. 

As a matter of comparison we recall that the "collide 
and stream" LB model is equivalent to an upwind FDLB 
model, when also the relaxation term is calculated on 
the characteristics line |31|] and the choice 6s — c5t is 
adopted. The resulting apparent value of the viscosity 



is i^ap = xc (t — St/2). For this reason the "collide and 
stream" LB model suffers mainly from the lack of sta- 
bility when T ~ (5^/2 so that very low values of viscosity 
cannot be accessed |25j. 



B. Flux limiter schemes 

Figure n shows two characteristics lines on the square 
lattice involving the distribution functions fi{x,t) and 
/^{Xjt), respectively. For convenience, we denote g^j the 
value of the quantity gi in node j at time t = kSt. Ac- 
cording to the general procedure to construct high order 
Total Variation Diminishing (TVD) schemes using flux 
limiters 0, we rewrite the updating rule lfT7|l in 

a conservative form using two fluxes [s^. l55l Is^ 



fk+l _ J-fc _ £^ 



+ 5tQi, (26) 



where 
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and 



■f^j,0-l) + l/2 



(28) 



The flux limiter '^{Qfj) introduced in (|27|) is expressed 
as a function of the smoothness 
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In particular, the second order Lax - Wendroff scheme 
is recovered for '^{Qf^) = 1. The upwind scheme, de- 
scribed in the previous subsection, is recovered as another 
particular case, when '^{Qfj) = 0. A wide choice of flux 
limiters are at our disposal in the literature 0, . 
LB simulations reported in this paper were done using 
the Monitorized Central Difference (MCD) hmiter |4J|] 
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but other limiters give qualitatively similar results. 

Eqs. 1)26(1 satisfy the global particle and momentum 
conservation. When using the flrst order upwind scheme, 
the spurious terms introduced in the mass and momen- 
tum conservation equations are linearly dependent on the 
lattice spacing ds. Since flux limiter schemes are adapt- 
ing themselves to the local smoothness (|29|) of the dis- 
tribution functions, it is rather cumbersome to derive 
analytical expressions of the spurious numerical term tp 
in these cases. LB simulations of diffusion phenomena 
done using flux limiter schemes suggest a second order 
dependence of the value V' on the lattice spacing 5s [s^ 
such that ip = {5sY /2cL and the apparent value of the 
kinematic viscosity (I25II becomes 



^ap_flux — 



2cL 



2 1 



(31) 



(27) 



When the lattice spacing is a small quantity, the use of 
flux limiter schemes is expected to improve the accuracy 
of FDLB simulations as well as the stability. 

A different approach that allows to avoid spurious 
terms in the conservation equations is provided by the 
volumetric LB scheme introduced in Ref. 46] which sat- 
isfies detailed balance and achieves the desired order of 
accuracy. We will refer to the fractional version of the 
aforementioned scheme constructed in the case of a ho- 
mogeneous fluid on a uniform mesh ^3 since we are using 
a regular and uniform lattice. In the scheme proposed in 
Ref. ^3 the value of the viscosity can be reduced with 
respect to the " collide and stream" LB and the Courant- 
Friedrichs-Levy number CFL = cSt/ds can be smaller 
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than 1. Moreover, some unphysical spurious invariants 
are removed. Our scheme can have very small values of 
viscosity since the numerical contribution to the value of 
viscosity, proportional to ip, can be reduced and made 
much smaller than the physical term, proportional to r, 
without having stability problems. This depends on the 
fact the in finite difference schemes the values of 5s and 
5t can be set independently from the value of c. Our 
choice of 5t and 5s is such to guarantee that the CFL 
number is much less than 1 and that the unavoidable 
spurious terms, introduced by the numerical scheme and 
proportional to (-0 — , can be done as small as desired. 



IV. SIMULATION RESULTS 

In this Section we report the results of our simula- 
tions. For all runs we used N = 1024, 5s = 1/256, 
and 5t — 10~^. In the following, lengths are expressed 
in units of lattice spacing and time is expressed as the 
product of algorithm steps by 5t. All quenches below the 
critical temperature were to the temperature 9 = 0.79 
where the coexisting densities are nuquid — 1.956 and 
TT-vapor = 0.226. Each simulation was started with small 
fluctuations (0.1%) in the density about a mean value 
n that was either symmetric {n = 1.09, liquid fraction 
f3 = 0.5) or slightly off-symmetric {n = 1.0, f3 = 0.45). 
The parameter k controlling the surface tension was set 
to 5 X 10~^ to have an interface thickness of ~ 6 lattice 
spacings. The viscosity was varied by changing t. We 
fixed an upper bound of r by the following argument. 
It is well known that the continuum hypothesis and the 
Navier - Stokes equation are valid only for small values 
of the Knudsen number Kn |58|. Since Kn — ct/L and 

c ~ L ~ 1 in our simulations, this means r ~ 10~^. We 
implemented the upwind and the flux limiter schemes 
and compared results when r = 10~^. In this case the 
spurious numerical contribution of the upwind scheme is 
larger than the physical one. Numerical contributions 
get negligible when the flux limiter scheme is considered 
instead. Therefore one expects to observe qualitative and 
quantitative differences. We used also the value r = 10^'^ 
with the flux limiter scheme to access a higher viscosity 
regime. 

In order to have different and independent tools to es- 
timate the domains size we used the following quanti- 
ties: Ri{t), the inverse of the length of the interfaces 
of domains, measured by counting lattice points where 
the order parameter p{x,t) — n(x,t) — n is such that 
p{x,t)p{x' ,t) < 0; R2{t), the inverse of the first moment 
of the spherically averaged structure factor 



R2{t) 



JC{k,t)dk 
JkC{k,t)dk' 



(32) 



where k = \k\is the modulus of the wave vector in Fourier 
space, and 
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FIG. 2: Evolution of domains size recovered for r = 10~ , 
/3 = 0.5 with the flux limiter scheme: Ri (A), R2 (o), 
(•). R's are measured in lattice spacings and Ri has been 
multiplied by 4,000,000 to be shown in the same plot. The 
straight line has slope 2/3. 



with p{k,t) the spatial Fourier transform of the order 
parameter p{x,t). The angle brackets denote an average 
over a shell in k space at fixed k. The last quantity 
i?3(t) is defined as the inverse of the first moment of the 
spherically averaged structure factor of the fluid velocity 
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(34) 
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(33) 



with Cu{k,t) —< \u{k)\'^ >. In all the figures Ri was 
multiplied by 4,000,000 to be shown in the same plot 
with i?2 and R3. 

In Fig. 121 we present the three measures of domains size 
as function of time for the case r = 10^^ with flux limiter 
scheme, and symmetric composition. It is interesting to 
note that R2 and R3 have the same trend, with a similar 
prefactor. This feature holds for all the runs we consid- 
ered. This last point is not obvious a priori. After a swift 
initial growth the evolution of all quantities suggests the 
existence of the growth exponent 2/3. This is in accor- 
dance with previous studies on symmetric liquid-vapor 
systems at low viscosity when hydrodynamic flow is 
operating. In this regime hydrodynamics is the mech- 
anism to get domains circular since the flow is driven 
by the difference in Laplace pressure between points of 
different curvature on the boundary of domains. This re- 
mark is confirmed when looking at configurations of the 
density n. In Fig. |3| we show contour plots of a part of 
the whole system at consecutive times. The vapor bub- 
ble in the down left corner at t = 12, while evaporating, 
is rounded by the flow as it can be seen by comparing it 
with the shape at t = 15. 
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t=12 t=15 



FIG. 3: Contour plots of a portion 512 x 512 of the whole 
lattice of the density n in the case withr = lO""*, P = 0.5 and 
flux limiter scheme. Color code: black/white ^ liquid/vapor. 




wavelength 



FIG. 4: Velocity structure factor Cu(k) at time f = 15 in the 
case with r = 10"*, /3 = 0.5 and flux limiter scheme. Cu{k) 
is in arbitrary units and the wavelength is measured in lattice 
spacings. 



An indication about the velocity field comes from the 
structure factor Cu{k,t). In Fig. ^ we plot it at time 
t — 15. It exhibits a structure at a scale comparable 
with system size. All velocity components decay becom- 
ing small at low wavelengths and contributing little to 
the overall dynamics. A small bump can be seen at wave- 
length ~ 8 corresponding to capillary motion at interface 
length scale. A similar behavior was observed in binary 

fluids [13. 

In the case with the upwind scheme the estimation of 
the growth exponent is more difficult since data are noisy 
and none of the i?'s shows a clear trend. From Fig. it 
seems the system enters a late regime characterized by 
an exponent consistent with the value 1/2. We believe 
that this behavior is due to spurious terms in the macro- 
scopic equations that are considerably larger when using 
the upwind scheme than in the case with flux limiter. 
These terms produce a numerical diffusivity when they 
are not negligible. This is confirmed by the analysis of 
the velocity fields in the two cases. In Fig. |Hlwe plot the 
order parameter p and velocity modulus u along a hori- 
zontal cross section of the system taken at the same long 
time. Two comments are in order here. It is quite un- 
avoidable to have spurious velocities at interfaces where 
density gradients are present with LB models (irrespec- 
tively of the particular model used H^l). And also the 
present model shows this unpleasant feature. Nonethe- 
less it is evident that the flux limiter scheme allows to 
dump considerably these spurious contributions. Indeed, 
with flux limiter the maximum value of velocity at inter- 
face is 0.13 {Ma = u/cs = 0.14) while with the upwind it 
is about 2 times larger being 0.23 (Ma = 0.26). The high 
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FIG. 5: Evolution of domains size in the case with t = 10~ , 
/3 = 0.5 and the upwind scheme: Ri (A), R2 (o), R3 (•). R's 
are measured in lattice spacings and Ri as been rescaled by 
4,000,000 to be shown in the same plot. The straight line has 
slope 1/2. 



value of the Mach number Ma makes the expansions (0) 
less reliable with the upwind scheme. 

Due to the better performance of the flux limiter 
scheme we decided to adopt it for further simulations. 
In Fig. [7| we plot the three measures of domains size as 
function of time for the case r = 10"'^ with symmet- 
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FIG. 6: Order parameter p (dashed line) and velocity modulus 
u (full line) profiles are shown along the line taken at y = 256 
lattice spacings from bottom at time t = 20 for the upwind 
scheme (upper panel) and flux limiter scheme (lower panel) 
with r = 10"'', f3 = 0.5. 



ric composition. After initial growth all the quantities 
suggest the existence of the growth exponent 1/2. This 
is in accordance with previous studies on liquid-vapor 
systems at high viscosity [4l| at symmetric composition 
when growth is expected to be described by the AUen- 
Cahn theory of interfaces dynamics which gives an expo- 
nent 1/2 |38| and hydrodynamics is not operating. Due 
to limits imposed by system size we cannot access very 
long times to probe whether the hydrodynamic regime is 
the late regime as previously argued |43| . Fig. shows 




FIG. 7: Evolution of domains size in the case with r = 10~ , 
/3 = 0.5 and flux limiter scheme: Ri (A), R2 (o), R3 (•). R'a 
are measured in lattice spacings and Ri as been multiplied 
by 4,000,000 to be shown in the same plot. The straight line 
has slope 1/2. 




t=12 t=15 



FIG. 8: Contour plots of a portion 512 x 512 of the whole 
lattice of the density n in the case with r — 10"^, (3 = 0.5 and 
flux limiter scheme. Color code: black/white liquid/vapor. 

density contour plots at consecutive times. Growth seems 
to be mainly driven by evaporation of vapor domains. 

Finally, we considered the case of an off-symmetric sys- 
tem with a liquid fraction f3 — 0.45. In Fig. we plot 
the three measures of domains size as function of time 
for the case r = 10~*. After the initial growth all the 





t 

FIG. 9: Evolution of domains size in the case with r — 10~^, 
l3 = 0.45 and flux limiter scheme: Ri (A), R2 (o), _R3 (•). 
R's are measured in lattice spacings and Ri as been multiplied 
by 4,000,000 to be shown in the same plot. The straight lines 
have slopes 2/3 and 1/2. 



quantities suggest the existence of a growth exponent 2/3 
which quite soon changes to 1/2. In previous studies 
with a free-energy LB model it was found the growth ex- 
ponent to be 2/3 with liquid fraction (3 = 0.31 and 
2/3 crossing over to 1/2 at /3 = 0.17 113. The prob- 
lem of off-symmetric liquid-vapor phase separation was 
recently addressed in Ref. (61|. There it was pointed out 
that in the case of a dispersion of liquid drops in vapor, 
the growth should proceed with an exponent 1/2 and the 
result was proven in the case of highly asymmetric com- 
position with P = 0.1. We believe that we are probing 
a regime similar to that seen in two-dimensional binary 
fluids where, once hydrodynamics flow has made domains 
circular, AUen-Cahn growth takes over 36] . This inter- 
pretation seems to be confirmed by configurations of the 
system, presented in Fig.^| They show that liquid drops 
in the vapor matrix are almost circular at i = 6 so that 
the hydrodynamic mechanism is no more effective. 



V. CONCLUSIONS 

The correct choice of the numerical scheme is essential 
to recover the real physics of a fluid system subjected 
to LB simulations. In the case of a liquid- vapor system 
we have seen that simulation results exhibit significant 
changes when the numerical contribution of the finite dif- 
ference scheme to the apparent value of the transport co- 
efficients becomes comparable with the expected physical 
value. The numerical contribution of the first order up- 
wind scheme is linearly dependent on the lattice spacing 



t=6 t=10 

FIG. 10: Contour plots of a portion 512 x 512 of the whole 
lattice of the density n in the case with r = 10"*, 13 = 0.45. 
Color code: black/white — > liquid/vapor. 



5s and switches to an higher order for the flux limiter 
scheme. Since Ss < 1, the flux limiter scheme reduces 
the computing effort in terms of required lattice nodes 
and gives physical results which are more accurate for 
the same number of lattice nodes per unit length. Spu- 
rious velocities at interfaces can be considerably damped 
and very low viscosity systems can be simulated preserv- 
ing numerical stability. The main limitation comes from 
the requirement of a small time step when very low val- 
ues of viscosity are needed. To give an idea of the CPU 
time we report that our code takes 6 hours to perform 
10^ algorithm steps by using 32 Xeon 3.055 GHz pro- 
cessors on the IBM Linux Cluster 1350 at CINECA |63| 
with Myrinet IPC network and the Portable Extensible 
Toolkit for Scientific Computation (PETSc 2.1.6) devel- 
OT^d at Argonne National Laboratory, Argonne, Illinois 

The model allowed to clarify the picture of phase sepa- 
ration in liquid-vapor system. We found that the growth 
exponent depends on either the fluid viscosity and the 
system composition. When liquid and vapor are present 
in the same amount, the growth exponent is 2/3 and 
1/2 at low and high viscosity, respectively. When the 
liquid fraction is less abundant than the vapor one, we 
can access a late time regime at low viscosity. In this 
regime the hydrodynamic transport is no more effective 
so we are able to see the crossover from the exponent 2/3 
to 1/2 which is characteristic of the AUen-Cahn growth 
mechanism. 

Finally, we note that our results as well as previous 
ones have been obtained in the case of isothermal sys- 
tems. It would be interesting to incorporate the energy 
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conservation equation into the model to allow non uni- 
form temperatures in the liquid- vapor system undergoing 
phase separation. 
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